Using the sat dataset, fit a model with the total SAT score as the response and expend, salary, ratio and takers as predictors. Perform regression diagnostics on this model to answer the following questions. Display any plots that are relevant. Do not provide any plots about which you have nothing to say. Suggest possible improvements or corrections to the model where appropriate.
n <- dim(sat)[1]
p <- dim(sat)[2]
model <-
lm(total ~ expend + salary + ratio + takers, data = sat)
sat %>%
select(total, expend, salary, ratio, takers) %>%
ggpairs(progress = F)
Lo primero que hago es crear una matriz de dispersión para tener una idea de como se distribuyen las variables y la relación 2x2 que tienen entre si. Utilizo la librería GGally que tiene la función ggpairs.
par(mfrow = c(2,2))
plot(model)
Gráficos Diagnósticos
Para evaluar la homocedasticidad del error utilizo los gráficos de:
En este caso se observa cierta irregularidad en el gráfico \(\hat{y} \times \hat{\epsilon}\) que podría traducir no-linearidad, aunque dentro de todo no se observa un patrón definido, ni tan anormal que haga pensar que el erro no sigue una homocedasticidad.
shapiro.test(resid(model)) %>%
pander()
| Test statistic | P value |
|---|---|
| 0.9769 | 0.4304 |
Para evaluar la normalidad del error, utilizo el gráfico de Q-Q de los residuales. En este caso el modelo se acerca bastante a la linea recta, con algunos valores raros que se alejan, que corresponden a los con valores mas altos de residuales, y que habría que examinar mejor.
Por otra parte podemos utilizar el test de Shapiro-Wilk cuya H0 es que el error no difiere de la distribución normal. En este caso no hay evidencia para rechazar la H0
halfnorm(hatvalues(model), labs = rownames(sat))
augment(model) %>%
filter(.hat > 2*(p/n)) %>%
pander()
| .rownames | total | expend | salary | ratio | takers | .fitted | .se.fit | .resid | .hat | .sigma | .cooksd | .std.resid |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| California | 902 | 4.992 | 41.08 | 24 | 45 | 917.8 | 17.37 | -15.85 | 0.2821 | 32.95 | 0.02571 | -0.572 |
| Utah | 1076 | 3.656 | 29.08 | 24.3 | 4 | 1010 | 17.67 | 65.77 | 0.2921 | 30.9 | 0.4715 | 2.39 |
Para evaluar los valores extremos en el modelo, utilizamos dos gráficos.
En este caso podemos ver que los valores para California y Utah cumplen con las características para valores extremos.
(outliers <-
augment(model) %>%
mutate(.t.resid = rstudent(model)) %>%
arrange(-abs(.t.resid)) %>%
head(5)) %>%
pander()
| .rownames | total | expend | salary | ratio | takers | .fitted | .se.fit | .resid | .hat | .sigma | .cooksd |
|---|---|---|---|---|---|---|---|---|---|---|---|
| West Virginia | 932 | 6.107 | 31.94 | 14.8 | 17 | 1023 | 8.147 | -90.53 | 0.06207 | 29.92 | 0.1081 |
| Utah | 1076 | 3.656 | 29.08 | 24.3 | 4 | 1010 | 17.67 | 65.77 | 0.2921 | 30.9 | 0.4715 |
| North Dakota | 1107 | 4.775 | 26.33 | 15.3 | 5 | 1040 | 9.31 | 66.57 | 0.08105 | 31.37 | 0.07954 |
| New Hampshire | 935 | 5.859 | 34.72 | 15.6 | 70 | 869.1 | 9.412 | 65.87 | 0.08284 | 31.4 | 0.07989 |
| Nevada | 917 | 5.16 | 34.84 | 18.7 | 30 | 971.1 | 6.967 | -54.15 | 0.04539 | 32 | 0.02731 |
| .std.resid | .t.resid |
|---|---|
| -2.859 | -3.124 |
| 2.39 | 2.53 |
| 2.124 | 2.214 |
| 2.103 | 2.19 |
| -1.695 | -1.732 |
abs(outliers$.t.resid) > abs(qt(.05/(n * 2), n - p))
## [1] FALSE FALSE FALSE FALSE FALSE
Para valorar posibles outliers podemos utilizar la estrategia planteada en el libro de Faraway, y calcular los residuales studentizados. En este sentido los residuales cuyo valor superen un limite determinado por la corrección de Bonferroni, se podrían considerar posibles outliers.
En este caso el valor que mas puede ser un outlier es West Virginia, aunque el valor de .t.resid de west virgina no es mayor al valor limite de p corregido.
model_1 <-
sat %>%
rownames_to_column() %>%
filter(rowname != "West Virginia") %>%
lm(total ~ expend + salary + ratio + takers, data = .)
pander(model)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 1046 | 52.87 | 19.78 | 7.858e-24 |
| expend | 4.463 | 10.55 | 0.4231 | 0.6742 |
| salary | 1.638 | 2.387 | 0.6861 | 0.4962 |
| ratio | -3.624 | 3.215 | -1.127 | 0.2657 |
| takers | -2.904 | 0.2313 | -12.56 | 2.607e-16 |
pander(model_1)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 1058 | 48.51 | 21.8 | 3.259e-25 |
| expend | 7.359 | 9.693 | 0.7592 | 0.4518 |
| salary | 1.088 | 2.191 | 0.4963 | 0.6221 |
| ratio | -3.94 | 2.943 | -1.338 | 0.1876 |
| takers | -2.972 | 0.2127 | -13.98 | 8.685e-18 |
Como West Virgina es un valor extremo (esta lejos de la linea del 0 en el plot de \(\hat{y} \times \sqrt{|{\hat{\epsilon}}|}\))
plot(model, 4)
halfnorm(cooks.distance(model), labs = rownames(sat))
model_2 <-
sat %>%
rownames_to_column() %>%
filter(rowname != "Utah") %>%
lm(total ~ expend + salary + ratio + takers, data = .)
pander(model)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 1046 | 52.87 | 19.78 | 7.858e-24 |
| expend | 4.463 | 10.55 | 0.4231 | 0.6742 |
| salary | 1.638 | 2.387 | 0.6861 | 0.4962 |
| ratio | -3.624 | 3.215 | -1.127 | 0.2657 |
| takers | -2.904 | 0.2313 | -12.56 | 2.607e-16 |
pander(model_1)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 1058 | 48.51 | 21.8 | 3.259e-25 |
| expend | 7.359 | 9.693 | 0.7592 | 0.4518 |
| salary | 1.088 | 2.191 | 0.4963 | 0.6221 |
| ratio | -3.94 | 2.943 | -1.338 | 0.1876 |
| takers | -2.972 | 0.2127 | -13.98 | 8.685e-18 |
pander(model_2)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 1094 | 53.42 | 20.48 | 4.043e-24 |
| expend | -0.9427 | 10.19 | -0.0925 | 0.9267 |
| salary | 3.096 | 2.328 | 1.33 | 0.1904 |
| ratio | -7.639 | 3.428 | -2.229 | 0.031 |
| takers | -2.931 | 0.2188 | -13.4 | 3.946e-17 |
Para evaluar los puntos mas influyentes podemos utilizar la distancia de Cook que refleja los cambios en el modelo al no incluir una observación determinada. En este caso destacan Utah y West Virginia como observaciones influyentes. Además podemos volver a mirar el gráfico de Residuales v/s leverage, que también refleja puntos con un grado alto de palanca y un residual alto.
influence(model)$coefficients %>%
as_tibble(rownames = "id") %>%
gather(key, value, -id) %>%
ggplot(aes(fct_rev(id), value, group = key)) +
geom_line(alpha = .75)+
coord_flip() +
geom_hline(yintercept = 0, linetype = 3, color = "darkred") +
facet_grid(~ key, scales = "free") +
scale_colour_viridis_d(end=.85)
La función influence permite obtener los coeficientes al restar una observación \(x_i\) (“leave one out”). De esta manera puedo construir un gráfico que revela los cambios en los diferentes coeficientes para cada una de las observaciones. En este gráfico se ve claramente como Utah genera los mayores cambios en todos los coeficientes.
plot(model, 5)
pander(sat[44,])
| expend | ratio | salary | takers | verbal | math | total | |
|---|---|---|---|---|---|---|---|
| Utah | 3.656 | 24.3 | 29.08 | 4 | 513 | 563 | 1076 |
El modelo multivariable es explicativo, con un valor de R^2 relativamente alto. En el modelo el gasto por estudiante expend y el sueldo salary se relacionan de manera positiva con el resultado del test total. La ratio de profesores y el porcentaje de estudiantes que hacen el examen de manera negativa.
Utah es un caso raro en la medida que tiene un score muy alta en relación a un gasto expend y un sueldo salary inferior al promedio por una parte, y una ratio. El porcentaje de alumnos elegibles takers llama la atención, ya que parece extremadamente bajo en relación a los otros parámetros.
En primer lugar me remito al gráfico de dispersión de (a), que permite reflejar la distribución y la relacion 1 a 1 de las variables.
# added variable
car::avPlots(model)
## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4
Para evaluar la estructura del modelo son de utilidad los gráficos de regresión parcial. La regresión parcial permite aislar la influencia de una variable independiente con la variable dependiente, restando la influencia de las otras variables independientes. Además permiten valorar la presencia de outliers y puntos influyentes.
En este caso podemos ver la clara relación inversa de la variable takers con el total y la relación mas débil de las otras variables. POdemos observar ademas el caso raro de Utah
par(mfrow = c(2,2))
termplot(model, partial.resid = T)
Los gráficos de residuales parciales, permiten conocer la respuesta aislando el efecto del resto de variables independientes. También permiten poder observar la presencia de diferencias estructurales, de no -linearidad
Using the teengamb dataset, fit a model with gamble as the response and the other variables as predictors. Answer the questions posed in the previous question.
n <- dim(teengamb)[1]
p <- dim(teengamb)[2]
teengamb_cl <-
teengamb %>%
rownames_to_column(".rownames") %>%
mutate(.rownames = factor(.rownames),
sex = factor(sex, levels = c (0,1), labels = c('male', 'female')))
model <-
lm(gamble ~ sex + status + income + verbal, data = teengamb_cl)
teengamb_cl %>%
select(gamble, sex, status, income, verbal) %>%
ggpairs(progress = F, aes(color = sex), lower = list(combo = wrap("facethist", bins = 10)))
En este caso gráfico con utilizando sex como un factor.
par(mfrow = c(2,2))
plot(model)
Gráficos Diagnósticos
En este caso el gráfico podría definir cierta heterocedasticidad en la medida que a mayor valor de \(\hat{y}\) parece haber mayor dispersión del error.
var.test(resid(model)[teengamb_cl$sex == "female"], resid(model)[teengamb_cl$sex == "male"]) %>%
pander()
| Test statistic | num df | denom df | P value | Alternative hypothesis | ratio of variances |
|---|---|---|---|---|---|
| 0.316 | 18 | 27 | 0.01374 * | two.sided | 0.316 |
En este caso al menos en cuanto a las poblaciones según sexo, hay una diferencia significativa en las varianzas
shapiro.test(resid(model)) %>%
pander()
| Test statistic | P value |
|---|---|
| 0.8684 | 8.16e-05 * * * |
En este caso podemos observar que el gráfico Q-Q difiere de la normalidad con un patrón de cola larga. Por otra parte la preba de Shapiro Wilk refleja que se rechaza la H0 de normalidad de la distribución del error.
halfnorm(hatvalues(model), labs = teengamb_cl$.rownames)
augment(model) %>%
mutate(.rownames = teengamb_cl$.rownames) %>%
filter(.hat > 2*(p/n)) %>%
pander()
| gamble | sex | status | income | verbal | .fitted | .se.fit | .resid | .hat | .sigma | .cooksd | .std.resid | .rownames |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 88 | male | 18 | 12 | 2 | 77.12 | 11.1 | 10.88 | 0.2395 | 22.88 | 0.01904 | 0.5498 | 31 |
| 90 | male | 38 | 15 | 7 | 78.25 | 10.68 | 11.75 | 0.2213 | 22.87 | 0.01957 | 0.5867 | 33 |
| 14.1 | male | 28 | 1.5 | 1 | 28.5 | 12.67 | -14.4 | 0.3118 | 22.8 | 0.05304 | -0.7651 | 35 |
| 69.7 | male | 61 | 15 | 9 | 73.54 | 12.46 | -3.836 | 0.3016 | 22.95 | 0.003535 | -0.2023 | 42 |
En este caso podemos ver que los valores para 42 y 35 cumplen con las características para valores extremos. Habría que revisar también los 31 y 33.
augment(model) %>%
mutate(.t.resid = rstudent(model)) %>%
filter(abs(.t.resid) > abs(qt(.05/(n * 2), n - p))) %>%
pander()
| gamble | sex | status | income | verbal | .fitted | .se.fit | .resid | .hat | .sigma | .cooksd | .std.resid | .t.resid |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 156 | male | 27 | 10 | 4 | 61.75 | 7.984 | 94.25 | 0.1238 | 16.74 | 0.5565 | 4.438 | 6.016 |
En este caso el valor 24 parece ser un outlier
model_1 <-
teengamb_cl %>%
filter(.rownames != "24") %>%
lm(gamble ~ sex + status + income + verbal, data = .) %>%
pander()
pander(model)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 22.56 | 17.2 | 1.312 | 0.1968 |
| sexfemale | -22.12 | 8.211 | -2.694 | 0.01011 |
| status | 0.05223 | 0.2811 | 0.1858 | 0.8535 |
| income | 4.962 | 1.025 | 4.839 | 1.792e-05 |
| verbal | -2.959 | 2.172 | -1.362 | 0.1803 |
pander(model_1)
## Warning in pander.default(model_1): No pander.method for "knit_asis",
## reverting to default.
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 7.631 | 12.93 | 0.5904 | 0.5582 |
| sexfemale | -16.3 | 6.133 | -2.657 | 0.01118 |
| status | 0.1739 | 0.2083 | 0.8346 | 0.4088 |
| income | 4.331 | 0.7636 | 5.672 | 1.265e-06 |
| verbal | -1.802 | 1.614 | -1.117 | 0.2707 |
La exclusión de la observación 24 genera varios cambios en todos los coeficientes del modelo. Ademas de ser un candidato a outlier es un valor muy influyente
plot(model, 4)
halfnorm(cooks.distance(model), labs = teengamb_cl$.rownames)
El gráfico de distancias de Cook confirma que el valor 24 se aleja mucho del comportamiento de las otras variables, que es un valor influyente y probablemente un outlier.
influence(model)$coefficients %>%
as_tibble(rownames = "id") %>%
mutate(id = factor(id, levels = c(1:n))) %>%
gather(key, value, -id) %>%
ggplot(aes(fct_rev(id), value, group = key)) +
geom_line(alpha = .75)+
coord_flip() +
geom_hline(yintercept = 0, linetype = 3, color = "darkred") +
facet_grid(~ key, scales = "free")
plot(model, 5)
pander(model)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 22.56 | 17.2 | 1.312 | 0.1968 |
| sexfemale | -22.12 | 8.211 | -2.694 | 0.01011 |
| status | 0.05223 | 0.2811 | 0.1858 | 0.8535 |
| income | 4.962 | 1.025 | 4.839 | 1.792e-05 |
| verbal | -2.959 | 2.172 | -1.362 | 0.1803 |
teengamb_cl[24,] %>% pander()
| .rownames | sex | status | income | verbal | gamble | |
|---|---|---|---|---|---|---|
| 24 | 24 | male | 27 | 10 | 4 | 156 |
En primer lugar me remito al gráfico de dispersión de (a), que permite reflejar la distribución y la relación 1 a 1 de las variables.
car::avPlots(model)
En este caso el modelo refleja la relación directa del sexo masculino y del income con gamble, y la relación inversa suave con el score verbal. Es posible apreciar claramente que la observación 24 es un caso raro que se sale de las líneas de predicción para todos los coeficientes.
par(mfrow = c(2,2))
termplot(model, partial.resid = T)
En este caso los parciales de residuales reflejan lo que se ve en el primer gráfico de dispersión y es que hay variables que no tiene una distribución normal, como income y , que podrían hacer pensar en hacer que transformaciones de las variables mejorarían la estructura del modelo
For the prostate data, fit a model with lpsa as the response and the other variables as predictors. Answer the questions posed in the first question.
n <- dim(prostate)[1]
p <- dim(prostate)[2]
prostate_cl <-
prostate %>%
mutate(svi = as.logical(svi))
model <-
lm(lpsa ~ lcavol + lweight + age + lbph + svi + lcp + gleason + pgg45, data = prostate_cl)
ggpairs(prostate_cl, progress = F, lower = list(continuous = wrap("points", alpha = 0.5, size = .4)))
Lo primero que hago es crear una matriz de dispersión para tener una idea de como se distribuyen las variables y la relación 2x2 que tienen entre si. Utilizo la librería GGally que tiene la función ggpairs.
par(mfrow = c(2,2))
plot(model)
Gráficos Diagnósticos
Para evaluar la homocedasticidad del error utilizo los gráficos de:
En este caso no se observa un patrón definido por lo que se asume homocedasticidad.
shapiro.test(resid(model)) %>% pander()
| Test statistic | P value |
|---|---|
| 0.9911 | 0.7721 |
La gráfica Q-Q parece bastante pareja y aproximada a la linea. El test de Shapiro Wilk confirma la normalidad del error
halfnorm(hatvalues(model), labs = rownames(prostate))
augment(model) %>%
mutate(.rownames = rownames(prostate)) %>%
filter(.hat > 2*(p/n)) %>%
pander()
| lpsa | lcavol | lweight | age | lbph | svi | lcp | gleason | pgg45 | .fitted | .se.fit | .resid | .hat |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2.008 | 0.1823 | 6.108 | 65 | 1.705 | FALSE | -1.386 | 6 | 0 | 2.875 | 0.4072 | -0.867 | 0.3305 |
| 2.158 | 1.423 | 3.657 | 73 | -0.5798 | FALSE | 1.658 | 8 | 15 | 1.925 | 0.3311 | 0.2323 | 0.2184 |
| 2.298 | 0.6206 | 3.142 | 60 | -1.386 | FALSE | -1.386 | 9 | 80 | 2.049 | 0.3478 | 0.2481 | 0.241 |
| 3.075 | 1.839 | 3.237 | 60 | 0.4383 | TRUE | 1.179 | 9 | 90 | 3.544 | 0.3098 | -0.4689 | 0.1912 |
| 4.13 | 2.533 | 3.678 | 61 | 1.348 | TRUE | -1.386 | 7 | 15 | 4.07 | 0.3241 | 0.0593 | 0.2092 |
| .sigma | .cooksd | .std.resid | .rownames |
|---|---|---|---|
| 0.7034 | 0.1227 | -1.496 | 32 |
| 0.7119 | 0.004271 | 0.3709 | 37 |
| 0.7118 | 0.005703 | 0.402 | 41 |
| 0.7103 | 0.01423 | -0.736 | 74 |
| 0.7124 | 0.0002605 | 0.09413 | 92 |
Para evaluar los valores extremos en el modelo, utilzamos dos gráficos.
En este caso podemos ver que los valores para el caso 32, 37 y el 41 cumplen con las características para valores extremos.
augment(model) %>%
mutate(.t.resid = rstudent(model)) %>%
filter(abs(.t.resid) > abs(qt(.05/(n * 2), n - p))) %>%
pander()
| lpsa | lcavol | lweight | age | lbph | svi | lcp | gleason | pgg45 | .fitted |
| .se.fit | .resid | .hat | .sigma | .cooksd | .std.resid | .t.resid |
Para valorar posibles outliers podemos utilizar la estrategia planteada en el libro de Faraway, y calcular los residuales studentizados. En este sentido los residuales cuyo valor superen un limite determinado por la corrección de Bonferroni, se podrían considerar posibles outliers.
En este caso no hay valores que cumplan con lo propuesto
plot(model, 4)
halfnorm(cooks.distance(model), labs = rownames(prostate))
model_1 <-
prostate_cl %>%
rownames_to_column() %>%
filter(rowname != "32") %>%
lm(formula = lpsa ~ lcavol + lweight + age + lbph + svi + lcp +
gleason + pgg45, data = .)
model_2 <-
prostate_cl %>%
rownames_to_column() %>%
filter(rowname != "47") %>%
lm(formula = lpsa ~ lcavol + lweight + age + lbph + svi + lcp +
gleason + pgg45, data = .)
pander(model)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 0.6693 | 1.296 | 0.5163 | 0.6069 |
| lcavol | 0.587 | 0.08792 | 6.677 | 2.111e-09 |
| lweight | 0.4545 | 0.17 | 2.673 | 0.008955 |
| age | -0.01964 | 0.01117 | -1.758 | 0.08229 |
| lbph | 0.1071 | 0.05845 | 1.832 | 0.0704 |
| sviTRUE | 0.7662 | 0.2443 | 3.136 | 0.002329 |
| lcp | -0.1055 | 0.09101 | -1.159 | 0.2496 |
| gleason | 0.04514 | 0.1575 | 0.2867 | 0.775 |
| pgg45 | 0.004525 | 0.004421 | 1.024 | 0.3089 |
pander(model_1)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 0.1719 | 1.329 | 0.1293 | 0.8974 |
| lcavol | 0.5653 | 0.08847 | 6.39 | 7.93e-09 |
| lweight | 0.6217 | 0.202 | 3.077 | 0.002793 |
| age | -0.02127 | 0.01115 | -1.908 | 0.05963 |
| lbph | 0.09559 | 0.05853 | 1.633 | 0.106 |
| sviTRUE | 0.7604 | 0.2426 | 3.135 | 0.002347 |
| lcp | -0.106 | 0.09036 | -1.173 | 0.244 |
| gleason | 0.05069 | 0.1564 | 0.3241 | 0.7466 |
| pgg45 | 0.004468 | 0.00439 | 1.018 | 0.3116 |
pander(model_2)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 0.0488 | 1.29 | 0.03783 | 0.9699 |
| lcavol | 0.5659 | 0.08615 | 6.568 | 3.574e-09 |
| lweight | 0.4577 | 0.1657 | 2.762 | 0.007009 |
| age | -0.01685 | 0.01095 | -1.538 | 0.1277 |
| lbph | 0.1167 | 0.05711 | 2.043 | 0.04409 |
| sviTRUE | 0.8265 | 0.2394 | 3.452 | 0.0008625 |
| lcp | -0.09798 | 0.08876 | -1.104 | 0.2727 |
| gleason | 0.1142 | 0.1562 | 0.7313 | 0.4665 |
| pgg45 | 0.004457 | 0.004309 | 1.034 | 0.3038 |
Los casos mas influyentes son el 32 y el 47. Si sacamos del modelo la observación 32, el coeficiente de lweight cambia en \(.15\). Al quitar la observación 47, el modelo se preserva mejor.
influence(model)$coefficients %>%
as_tibble(rownames = "id") %>%
mutate(id = factor(id, levels = c(1:n))) %>%
gather(key, value, -id) %>%
ggplot(aes(fct_rev(id), value, group = key)) +
geom_line(alpha = .75)+
coord_flip() +
geom_hline(yintercept = 0, linetype = 3, color = "darkred") +
facet_grid(~ key, scales = "free")
En este gráfico se ve como 32 genera un cambio importante en lweight. Llama la atención ademas la observación 69 que genera cambios en varios coeficientes y que tiene una distancia de cook importante.
plot(model, 5)
prostate_cl[c(69,32, 47),] %>% pander()
| lcavol | lweight | age | lbph | svi | lcp | gleason | pgg45 | lpsa | |
|---|---|---|---|---|---|---|---|---|---|
| 69 | -0.4463 | 4.409 | 69 | -1.386 | FALSE | -1.386 | 6 | 0 | 2.963 |
| 32 | 0.1823 | 6.108 | 65 | 1.705 | FALSE | -1.386 | 6 | 0 | 2.008 |
| 47 | 2.728 | 3.995 | 79 | 1.879 | TRUE | 2.657 | 9 | 100 | 2.569 |
Los casos 69, 32 y 47 son casos influyente. El caso 32 ademas tiene mucha palanca. Pero no tengo mucha evidencia como para decir que son outliers.
car::avPlots(model)
La regresión parcial confirma que 32 y mas aún el69 son casos raros que se comportan de manera diferente a la mayoría en cuanto a las lineas de regresión parcial.
par(mfrow = c(2,2))
termplot(model, partial.resid = T)
Aquí destaca el caso 32 con un lweight muy alejado del resto, y la estrucutra de lbph (y en menor grado, de lcp, pgg45 y gleason) que no es para nada homogenea, lo que ya se veía en el grafico de dispersion.
For the swiss data, fit a model with Fertility as the response and the other variables as predictors. Answer the questions posed in the first question.
n <- dim(swiss)[1]
p <- dim(swiss)[2]
model <-
lm(Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality, data = swiss)
swiss %>%
ggpairs(progress = F, lower = list(continuous = wrap("points", alpha = 0.75, size = .8)))
Lo primero que hago es crear una matriz de dispersión para tener una idea de como se distribuyen las variables y la relación 2x2 que tienen entre si. Utilizo la librería GGally que tiene la función ggpairs.
par(mfrow = c(2,2))
plot(model)
Gráficos Diagnósticos
No se observa un patrón definido, por lo que se acepta una homocedasticidad del error.
shapiro.test(resid(model)) %>% pander()
| Test statistic | P value |
|---|---|
| 0.9889 | 0.9318 |
El gráfico *Q-Q se adapta bien a la linea. La prueba no permite rechazar la H0 de normalidad del error.
halfnorm(hatvalues(model), labs = rownames(sat))
augment(model) %>%
filter(.hat > 2*(p/n)) %>%
pander()
| .rownames | Fertility | Agriculture | Examination | Education | Catholic | Infant.Mortality | .fitted | .se.fit |
|---|---|---|---|---|---|---|---|---|
| La Vallee | 54.3 | 15.2 | 31 | 20 | 2.15 | 10.8 | 50.74 | 4.246 |
| V. De Geneve | 35 | 1.2 | 37 | 53 | 42.34 | 18 | 34.8 | 4.838 |
| .resid | .hat | .sigma | .cooksd | .std.resid |
|---|---|---|---|---|
| 3.562 | 0.3512 | 7.221 | 0.03437 | 0.6172 |
| 0.2024 | 0.4558 | 7.254 | 0.0002047 | 0.03829 |
Para evaluar los valores extremos en el modelo, utilizamos dos gráficos.
En este caso podemos ver que los valores para Vermont y Maine cumplen con las características para valores extremos. Llaman la atención también los puntos de La Valle y V. De Geneve que presentan valores extremos en una de las variables.
augment(model) %>%
mutate(.t.resid = rstudent(model)) %>%
filter(abs(.t.resid) > abs(qt(.05/(n * 2), n - p))) %>%
pander()
| .rownames | Fertility | Agriculture | Examination | Education | Catholic |
| Infant.Mortality | .fitted | .se.fit | .resid | .hat | .sigma | .cooksd | .std.resid |
| .t.resid |
No hay outliers por este criterio. Tampoco se observan en el gráfico de las distancias de Cook en los gráficos diagnósticos
par(mfrow = c(2,2))
plot(model, 4)
halfnorm(cooks.distance(model), labs = rownames(swiss))
influence(model)$coefficients %>%
as_tibble(rownames = "id") %>%
gather(key, value, -id) %>%
ggplot(aes(fct_rev(id), value, group = key)) +
geom_line(alpha = .75)+
coord_flip() +
geom_hline(yintercept = 0, linetype = 3, color = "darkred") +
facet_grid(~ key, scales = "free") +
scale_colour_viridis_d(end=.85)
colMeans(swiss) %>% pander()
| Fertility | Agriculture | Examination | Education | Catholic | Infant.Mortality |
|---|---|---|---|---|---|
| 70.14 | 50.66 | 16.49 | 10.98 | 41.14 | 19.94 |
swiss %>%
rownames_to_column() %>%
filter(rowname %in% c("Porrentruy", "Sierre")) %>%
pander()
| rowname | Fertility | Agriculture | Examination | Education | Catholic | Infant.Mortality |
|---|---|---|---|---|---|---|
| Porrentruy | 76.1 | 35.3 | 9 | 7 | 90.57 | 26.6 |
| Sierre | 92.2 | 84.6 | 3 | 3 | 99.46 | 16.3 |
Hay varias observaciones que resultan influyentes. POr le criterio de distancias de Cook, destaca Porrentruy que en segundo gráfico podemos ver, modifica agriculture y Infant.Mortality. Sierre modifica Infant.Mortality hacia el otro sentido, tiene una mortalidad baja, y Examination baja.
Hay dos puntos con alta palanca pero baja influencia, La Valle y V. De Geneve.
En primer lugar me remito al gráfico de dispersión de (a), que permite reflejar la distribución y la relación 1 a 1 de las variables.
car::avPlots(model)
La estructura parece correcta. Podemos ver como Sierre se aleja de la tendencia de los modelos.
par(mfrow = c(2,2))
termplot(model, partial.resid = T, se = T)
Se pueden observar los valores extremos pero con poca influencia en Education y en Infant.Mortality Destaca la estructura bimodal de Catholic. Se podría realizar una transformación de la variable a una categorial o incluso binaria.
Using the cheddar data, fit a model with taste as the response and the other three variables as predictors. Answer the questions posed in the first question.
n <- dim(cheddar)[1]
p <- dim(cheddar)[2]
model <-
lm(taste ~ Acetic + H2S + Lactic, data = cheddar)
cheddar %>%
ggpairs(progress = F)
par(mfrow = c(2,2))
plot(model)
Gráficos Diagnósticos
No se ve un patrón el el gráfico de valores ajustados v/s residuales que haga pensar en problemas de varianza del error.
shapiro.test(resid(model)) %>% pander()
| Test statistic | P value |
|---|---|
| 0.9802 | 0.8312 |
El gráfico Q-Q sigue una distribución muy cercana a la recta y la prueba de S-W no permite desechar la H0 por lo que se puede asumir la normalidad de los residuales.
halfnorm(hatvalues(model), labs = rownames(cheddar))
augment(model) %>%
rowid_to_column() %>%
filter(.hat > 2*(p/n)) %>%
pander()
| rowid | taste | Acetic | H2S | Lactic | .fitted | .se.fit | .resid | .hat |
| .sigma | .cooksd | .std.resid |
Los casos 20 y 26 son los casos mas extremos.
augment(model) %>%
mutate(.t.resid = rstudent(model)) %>%
filter(abs(.t.resid) > abs(qt(.05/(n * 2), n - p))) %>%
pander()
| taste | Acetic | H2S | Lactic | .fitted | .se.fit | .resid | .hat | .sigma |
| .cooksd | .std.resid | .t.resid |
No hay outliers según este criterio.
plot(model, 4)
halfnorm(cooks.distance(model), labs = rownames(cheddar))
influence(model)$coefficients %>%
as_tibble(rownames = "id") %>%
mutate(id = factor(id, levels = c(1:n))) %>%
gather(key, value, -id) %>%
ggplot(aes(fct_rev(id), value, group = key)) +
geom_line(alpha = .75)+
coord_flip() +
geom_hline(yintercept = 0, linetype = 3, color = "darkred") +
facet_grid(~ key, scales = "free")
Los casos 12 y 15 son casos influyentes.
model_1 <-
cheddar %>%
rownames_to_column() %>%
filter(rowname != "15") %>%
lm(formula = taste ~ Acetic + H2S + Lactic, data = .)
pander(model)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -28.88 | 19.74 | -1.463 | 0.1554 |
| Acetic | 0.3277 | 4.46 | 0.07349 | 0.942 |
| H2S | 3.912 | 1.248 | 3.133 | 0.004247 |
| Lactic | 19.67 | 8.629 | 2.28 | 0.03108 |
pander(model_1)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -17.76 | 17.62 | -1.007 | 0.3233 |
| Acetic | -2.47 | 4.004 | -0.617 | 0.5428 |
| H2S | 4.039 | 1.091 | 3.702 | 0.001061 |
| Lactic | 21.46 | 7.559 | 2.839 | 0.008862 |
El caso 15 modifica el coeficiente de Acetic en 2 , y Lactic en casi 2.
car::avPlots(model)
par(mfrow = c(2,2))
termplot(model, partial.resid = T)
NO se ven mayores problemas estructurales en los gráficos.
Using thehappydata, fit a model withhappyas the response and the other four variables aspredictors. Answer the questions posed in the first question.
Using thetvdoctordata, fit a model withlifeas the response and the other two variables aspredictors. Answer the questions posed in the first question.
For thedivusadata, fit a model withdivorceas the response and the other variables, except yearas predictors. Check for serial correlation.
Using the divusa data:
divorce as the response and unemployed, femlab, marriage, birth and military as predictors. Compute the condition numbers and interpret their meanings.# Miramos colinearidad
divusa %>%
select(-year, -divorce) %>%
cor() %>%
corrplot.mixed(., lower = "number", upper = "circle", tl.pos = "lt")
model <- lm(divorce ~ unemployed + femlab + marriage + birth + military, data = divusa)
Podemos observar que si hay varios predictores que presenta correlación entre si.
X <- model.matrix(model)[,-1]
eig <- eigen(t(X) %*% X)
eig$values
## [1] 1174600.548 21261.741 16133.842 6206.181 1856.894
(c_nums <- sqrt(eig$values[1] / eig$values))
## [1] 1.000000 7.432684 8.532498 13.757290 25.150782
Los números de condición son relativamente pequeños (< 30) por lo que a pesar de la correlación de parejas encontrada, no parece existir un problema grave de colinearidad.
VIF <- vif(X)
pander(rbind(vifs, SE = sqrt(vifs)))
Quitting from lines 789-791 (keymer_alejandro.Rmd) Error in rbind(vifs, SE = sqrt(vifs)) : object ‘vifs’ not found Calls:
Si bien hay algunas \(X_i\) donde la VIF se aleja de 1 que es la ortogonalidad, tampoco se aleja demasiado. Donde más se aleja es en \(X_{femlab}\) donde se puede interpretar que el error estándar aumenta en 1.9008618
model_1 <-
lm(formula = divorce+2*rnorm(dim(divusa)[1]) ~ unemployed + femlab + marriage + birth +
military, data = divusa)
pander(model)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 2.488 | 3.394 | 0.7331 | 0.4659 |
| unemployed | -0.1113 | 0.05592 | -1.989 | 0.05052 |
| femlab | 0.3836 | 0.03059 | 12.54 | 1.106e-19 |
| marriage | 0.1187 | 0.02441 | 4.861 | 6.772e-06 |
| birth | -0.13 | 0.0156 | -8.333 | 4.027e-12 |
| military | -0.02673 | 0.01425 | -1.876 | 0.06471 |
pander(model_1)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 1.956 | 5.229 | 0.374 | 0.7095 |
| unemployed | -0.06709 | 0.08616 | -0.7786 | 0.4388 |
| femlab | 0.3795 | 0.04713 | 8.054 | 1.33e-11 |
| marriage | 0.1229 | 0.03762 | 3.267 | 0.001674 |
| birth | -0.1327 | 0.02403 | -5.525 | 5.118e-07 |
| military | -6.257e-06 | 0.02195 | -0.000285 | 0.9998 |
A pesar de eso, si agregamos un poco de ruido el modelo se mantiene relativamente estable.
model_1 <-
lm(formula = divorce ~ unemployed + marriage + birth + military, data = divusa)
model_2 <-
lm(formula = divorce ~ unemployed + femlab + military, data = divusa)
# Modelo original
pander(model)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 2.488 | 3.394 | 0.7331 | 0.4659 |
| unemployed | -0.1113 | 0.05592 | -1.989 | 0.05052 |
| femlab | 0.3836 | 0.03059 | 12.54 | 1.106e-19 |
| marriage | 0.1187 | 0.02441 | 4.861 | 6.772e-06 |
| birth | -0.13 | 0.0156 | -8.333 | 4.027e-12 |
| military | -0.02673 | 0.01425 | -1.876 | 0.06471 |
# Modelo 1, divorce ~ unemployed + marriage + birth + military
pander(model_1)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 42.01 | 2.246 | 18.7 | 2.369e-29 |
| unemployed | -0.5686 | 0.07551 | -7.53 | 1.151e-10 |
| marriage | -0.05652 | 0.03566 | -1.585 | 0.1173 |
| birth | -0.2289 | 0.02396 | -9.552 | 1.978e-14 |
| military | -0.0155 | 0.02532 | -0.6121 | 0.5424 |
# VIF model 1
pander(vif(model_1))
| unemployed | marriage | birth | military |
|---|---|---|---|
| 1.295 | 1.927 | 1.925 | 1.245 |
# Modelo 2, divorce ~ unemployed + femlab + military
pander(model_2)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -3.784 | 1.232 | -3.072 | 0.002988 |
| unemployed | 0.01835 | 0.06086 | 0.3015 | 0.7639 |
| femlab | 0.4412 | 0.02412 | 18.3 | 5.529e-29 |
| military | -0.008315 | 0.02014 | -0.4128 | 0.681 |
# VIF model 2
pander(vif(model_2))
| unemployed | femlab | military |
|---|---|---|
| 1.276 | 1.075 | 1.195 |
En este caso reviso dos modelos alternativos. Uno sin la variable femlab cuyo VIF era alta en l modelo original, y un segundo modelo con la variable femlab pero sin marriage ni birth que tienen una correlación de pareja alta con femlab y entre si.
Los dos modelos alternativos pierden bastante poder explicativo lo que se ve reflejado en un \(R^2\) mas pequeño, sobretodo el sin femlab. Por otra parte, es verdad que en el primer modelo sin femlab no baja la VIF (aumenta en las variables que tiene correlación, birth y marriage). En el segundo modelo, sin birth y marriage, la \(R^2\) baja, pero la VIF también.
Como conclusión, creo que el modelo original es correcto.
For the longley data, fit a model with Employed as the response and the other variables as predictors.
model <- lm(Employed ~ ., data = longley)
X <- model.matrix(model)[,-1]
eig <- eigen(t(X) %*% X)
eig$values
## [1] 6.665299e+07 2.090730e+05 1.053550e+05 1.803976e+04 2.455730e+01
## [6] 2.015117e+00
(c_nums <- sqrt(eig$values[1] / eig$values))
## [1] 1.00000 17.85504 25.15256 60.78472 1647.47771 5751.21560
Los numeros de condición son bastante altos al igual que los eigenvalues, lo que sugiere multicolinearidad
corrplot.mixed(cor(X), lower = "number", upper = "circle", tl.pos = "lt")
En este caso las variables predictoras tienen mucha mayor correlación entre si.
VIF <- vif(X)
pander(rbind(vifs, SE = sqrt(vifs)))
Quitting from lines 872-874 (keymer_alejandro.Rmd) Error in rbind(vifs, SE = sqrt(vifs)) : object ‘vifs’ not found Calls:
En este caso excepto por Armed.Forces, y algo menos Unemployed, la mayoría de las variables predictoras correlacionan entre si, y hacen que aumente el error estándar del modelo. Es probable que un modelo mas simple pueda funcionar mejor.
pander(model)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -3482 | 890.4 | -3.911 | 0.00356 |
| GNP.deflator | 0.01506 | 0.08491 | 0.1774 | 0.8631 |
| GNP | -0.03582 | 0.03349 | -1.07 | 0.3127 |
| Unemployed | -0.0202 | 0.004884 | -4.136 | 0.002535 |
| Armed.Forces | -0.01033 | 0.002143 | -4.822 | 0.0009444 |
| Population | -0.0511 | 0.2261 | -0.2261 | 0.8262 |
| Year | 1.829 | 0.4555 | 4.016 | 0.003037 |
For the prostate data, fit a model with lpsa as the response and the other variables as predictors.
model <-
lm(lpsa ~ ., data = prostate)
X <- model.matrix(model)[,-1]
eig <- eigen(t(X) %*% X)
eig$values
## [1] 4.790826e+05 6.190704e+04 2.109042e+02 1.756329e+02 6.479853e+01
## [6] 4.452379e+01 2.023914e+01 8.093145e+00
(c_nums <- sqrt(eig$values[1] / eig$values))
## [1] 1.00000 2.78186 47.66094 52.22787 85.98499 103.73114 153.85414
## [8] 243.30248
HAlgunos valores son altos y otros no, lo que sugiere que podría haber multicolinerarida en varias combinaciones lineales.
corrplot.mixed(cor(X), lower = "number", upper = "circle", tl.pos = "d")
Se corrobora lo anterior en la medida de que hay alguno de los predictores que correlacionan bastante, como el caso de pgg45 con gleason y lcp; lcp con svl y lcavol.
VIF <- vif(X)
pander(rbind(vifs, SE = sqrt(vifs)))
Quitting from lines 919-921 (keymer_alejandro.Rmd) Error in rbind(vifs, SE = sqrt(vifs)) : object ‘vifs’ not found Calls:
De todas formas, a pesar de la correlación, no parece haber un mayor problema al analizar la VIF. La mayoría de los predictores suman poco error
Use the fat data, fitting the model described in Section 4.2.
Realizar el análisis completo de los residuos del modelo de regresión parabólico propuesto en lasección 1.2 con los datos de tráfico.2.
Realizar el análisis completo de los residuos de los modelos de regresión simple y parabólico pro-puestos en la sección 1.2 con los datos de tráfico, pero tomando como variable respuesta la velocidad(sin raíz cuadrada). Este análisis debe justificar la utilización de la raíz cuadrada de la velocidadcomo variable dependiente.